Universal reference state in a driven homogeneous granular gas 
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We study the dynamics of a homogeneous granular gas heated by a stochastic thermostat, in 
the low density limit. It is found that, before reaching the stationary regime, the system quickly 
"forgets" the initial condition and then evolves through a universal state that does not only depend 
on the dimensionless velocity, but also on the instantaneous temperature, suitably renormalized by 
its steady state value. We find excellent agreement between the theoretical predictions at Boltzmann 
equation level for the one-particle distribution function, and Direct Monte Carlo simulations. We 
conclude that at variance with the homogeneous cooling phenomenology, the velocity statistics 
should not be envisioned as a single-parameter, but as a two-parameter scaling form, keeping track 
of the distance to stationarity. 



I. INTRODUCTION 

A granular gas may be viewed as a collection of macroscopicparticles undergoing dissipative collisions. This very 
ingredient -inelasticity- gives rise to a rich phenomenology [J-Q, the understanding of which requires statistical 
mechanics tools: Kinetic theory has proven to be powerful at a microscopic/mesoscopic level of description, while 
at a macroscopic scale, hydrodynamic equations have been derived Q and put to the test. Yet, the relevance and 
consistency of a hydrodynamic framework, which is a central question, is still elusive 0-0] ■ 

Due to collisional dissipation, the granular temperature, defined as the variance of velocity fluctuations, decays 
monotonically in time in an isolated granular system Q. In the fast-flow regime, it has been shown numerically that 
for a wide class of initial conditions, the system reaches a homogeneous state in which all the time dependence of 
the one-particle distribution function is encoded in the temperature. This is the so-called homogeneous cooling state 
(HCS) which has been widely studied in the literature (9l[l0|. In such a homogeneous situation, the dynamics involves 
two time-scales: the kinetic or fast scale -a few collisions per particle- in which the scaling regime has not been 
reached yet and where the "microscopic" excitations relax, and the following hydrodynamic or slow scale in which the 
memory of initial conditions has been lost, and the velocity distribution evolves through the granular temperature 
[111 ]. Considering non homogeneous states, this separation of time scales opens the possibility of a hydrodynamic or 
coarse-grained description in terms of the density, velocity and tem per ature field, and the HCS then plays the role of 
the reference state when the Chapmann-Enskog method is applied [12j . 

On the other hand, several studies and experiments in granular matter deal with stationary states, which are reached 
under the action of some energy driving, often realized by a moving boundary, or by an interstitial medium that acts 
as a thermostat, see e.g. |13l - [l5| . In all these cases, the energy injected compensates for the energy lost in collisions. 
From a theoretical point of view, a minimalistic approach is to consider the system as driven by some random 
energy source, which can be implemented in different ways (l6j . For the hard particle model, one of the most used 
homogeneous heating is the so-called stochastic thermostat, which consists in a white noise force acting on each grain 
[lol [l7l - l27| . In the low density limit, the distribution function of the homogeneous state has been characterized [Ioj |. 
Hydrodynamic equations have been derived via the Chapmann-Enskog expansion (2lj and fluctuating hydrodynamics 
have been put forward in order to understand the large scale structure found in the stationary state [18| , or fluctuations 
of global quantities [25|. We stress that in all the studies pertaining to the hydrodynamics of a system heated by 
the stochastic thermostat [IH, [2l[, the stationary state played the role of "reference" state, as the HCS happens to 
be in the undriven case. It was therefore assumed that a one parameter scaling holds for the velocity probability 
distribution. The objective in this work is to analyze this point critically in the homogeneous case at the level of the 
Boltzmann equation. We will study, for arbitrary initial conditions, the type of state the system evolves into in a 
kinetic time scale. Surprisingly, we find that a universal state is reached in a kinetic scale -universal in the sense that 
it is independent of the initial conditions- but that depends on the quotient between the instantaneous temperature 
and its stationary value. 

The outline is as follows. Section [TT] opens with a definition of the model, and a summary of relevant previously 
known results. The key question addressed lies in the scaling form of the velocity distribution close to the steady 
state. Does it depend only on the suitably reduced velocity variable, as is the case in the HCS, or is another parameter 
relevant, that would encode the distance to stationarity ? We will argue in section Hi Al that a single parameter scaling 
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form is inconsistent. We shall then show in sections III Bl and III CI that a consistent two parameter scaling form can 
be identified. Its properties will be characterized by complementary numerical and analytical tools. Conclusions and 
perspectives will finally be discussed in section llTTl 



II. SCALING FORM OF THE VELOCITY DISTRIBUTION FUNCTION 



The system of interest is a dilute gas of N smooth inelastic hard particles of mass m and diameter a, which collide 
inelastically with a coefficient of normal restitution a independent of the relative velocity Jj|. If at time t there 
is a binary encounter between particles i and j, having velocities Vi(i) and Vj(i) respectively, the post-collisional 
velocities V ■ (t) and (t) are 



Vj = V i -l±^(*-V < j)*, 



2 

1 + a 



(<X • Vy)<T, 



(1) 



where V,j = V, — Vj is the relative velocity and & is the unit vector pointing from the center of particle j to the center 
of particle i at contact. Between collisions, the system is heated uniformly by a white noise acting independently on 



each grain [lfj, |18|, |20|-|22 



Fokkcr-Planck equation 



25| so that the one-particle velocity distribution, /(r, v, t), then obeys the Boltzmann- 
28f T For a homogeneous system this equation reads 



dv 2 ro(vi,v a )/(vi,t)/(v a ,t) + k^-f(y 1}t ), 



(2) 



where d is the dimension of space, £o measures the noise strength, and To is the binary collision operator 

r (vi, v a ) = / d<76(v 12 • 6-)(v 12 • <r)(a- 2 6- x - 1). (3) 

Here we have introduced the operator &" 1 which replaces the velocities Vi and v 2 by the precollisional ones and 
v 2 given by 



Vi = Vi 

v 2 = v 2 



1 + a 
2a 

1 + Q 

2a 



(& ■ Vi 2 )<T, 
(a ■ vi 2 )<r. 



(4) 



A. One-parameter scaling or beyond ? 

It is an observation from numerical simulations, that for a wide class of initial conditions the system reaches a 
stationary state [hJ HH, Hq[ ■ Assuming that total momentum is zero, i.e. J dw/(v,0) = 0, the state is characterized 
by an isotropic stationary distribution, f s (v). Let us define the scaled distribution function, Xs, by 

fs(v) = ^dXs(c), c=— , (5) 



where n is the density, v s = \j2T s jm is the thermal velocity and T s is the stationary temperature, tj 71 ^ = 
J dv^mv 2 f s (v). As Xs is rather close to a Maxwellian distribution, a reasonable strategy is to perform an ex- 
pansion in Sonine polynomials [29j . In the so-called first Sonine approximation, the steady state function then reads 

M 

Xs(c)~XM(c)[l+a^S 2 (c 2 )}, (6) 

where xm is the Maxwellian distribution with unit temperature, S^tc 2 ) = — ^pc 2 + ^c 4 , is the second Sonine 

polynomial, and af is the kurtosis of the distribution. Within this approximation, the distribution function can be 
calculated, with the result [l(| 

s( 16(l- Q )(l-2a 2 ) 
a2(a) 73 + 56d - 24da - 105a + 30(1 - a)a 2 ' lj 



3 



and a stationary temperature 



dv{d/2)e 



-i 2/3 



2ir 2 (1 — a 2 )na 



d-l 



(8) 



Now, let us consider an initial condition with a temperature that differs appreciably from the stationary temperature 
(we also assume that total momentum is zero, its precise value being immaterial). It is clear that the system will 
reach the stationary state in a hydrodynamic scale. The ensuing question is two-pronged. First of all, is the dynamics 
compatible with a universal scaling form -once the memory of initial condition is lost- that would provide a consistent 
solution to the Boltzmann equation, or is memory only washed out strictly speaking at the steady state point? Second, 
assuming such a scaling regime exists in some vicinity of the steady state, what is the minimal number of parameters 
required for its description? By analogy with unforced (HCS) phenomenology, a single parameter scaling might be 
anticipated: 



/(v,t) 



Xs(c) where c = v/vo(t) and Vo(t) = y / 2T(t)/m 



(9) 



is defined from the instantaneous temperature dnT(t) = J dvmv 2 /(v, t). We note that this scaling property holds for 
the Gaussian thermostat as well [l9| , where the particles are accelerated between collisions by a force proportional to its 
own velocity [3(j. Moreover, in the stochastic thermostat case, the one parameter scaling J[5| was implicitly assumed, 
and it seemed to be yield reasonable predictions at least close to the stationary state, see 0, HJ. Nevertheless, when 
the form ([9]) is inserted in the Boltzmann equation, Eq. ([2]), we obtain 



£o 2 d 2 
2«§(t) dej 



X.(ci) + 



1 dv {t) d 
Uo(f) dt dci 



[ciXs(ci)] 



d-l 



dc 2 To(ci,c 2 )x s (ci)x«(c2), 



(10) 



which is inconsistent: The left hand side depends on time while the right hand side does not. We conclude that such a 
solution should be ruled out, except when stationarity is reached and vq = v s . The hope is to capture the post- kinetic 
time dependence of /(v,i), through a more involved functional form, that would be free of the above inconsistency. 

As is customary, we shall seek for a normal solution J2f|, which in the present homogeneous case means that 
/(v, t) should only depend on time via the instantaneous granular temperature T(t). In conjunction with dimensional 
analysis, this leads to a function that should only depend on c and T(t)/T s = Vq/v 2 , which we write as 



/(v,t) 



X (c,/3), with (5 



v (t) 
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and again c = v/vo(t). Note that we have assumed isotropy, x depending on c = |c| and not on the full vector c, 
but this assumption can be easily relaxed. Note also that equivalent expressions can of course be chosen, such as 
X(v/v s ,(3). 

If a state such as (fTTT) holds, it represent a strong constraint on the form of the velocity distribution. The cor- 
responding dynamics can be partitioned in a first rapid stage -that we do not attempt to describe- where initial 
conditions matter, and a subsequent universal relaxation towards stationarity, where only the distance to the steady 
state is relevant, through the dimensionless inverse typical velocity ft — v s /vo(t). 



B. Numerical simulations answer 

To put the above scenario to the test, we have performed Direct Monte Carlo Simulations (DSMC) [3l| of N = 1000 
hard disks (d = 2) of unit mass and unit diameter, that collide inelastically with the collision rule given by equation 
([T]) . The thermostat is implemented following previous investigations fl8j and the results have been averaged over 10 5 
trajectories. For a given value of the inelasticity, we thus solve the time-dependent Boltzmann equation for different 
initial conditions and analyze whether, after some kinetic transient, all the time dependence of the distribution function 
goes through the dimensionless parameter ft. As it is difficult to measure the complete distribution function with 
the desired accuracy, we have worked with the cumulants of the scaled distribution, x( c j/?)- I n terms of the velocity 
moments, (v l ) = — J dvv l f(v,t), we have measured the kurtosis of the distribution 



(12) 
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which is proportional to the fourth cumulant of %(c, f$) and the quantity 



a 3 



(v 6 ) 3d {v A ) 



(d + 2)(d + 4) <u 2 ) 3 d+2{v 2 )' 



-2, 



(13) 



which can be viewed as the reduced sixth cumulant. If our scaling is correct, we expect that the cumulants quickly 
collapse for different initial conditions, as a function of (3. In Fig. [TJ we have plotted a-i and 03 versus /3 for a = 0.95. 
The initial conditions are either Maxwellian distributions with three different temperatures To, significantly above the 
steady state value T s , or asymmetric distributions made up of three possible velocities with different probabilities 



f(v x , v v , t = 0) = ^6 (v x + 8D/3) 8 (v v + 8D/3) + ^5 (v x - AD/3) S (v y - AD/ 3) 
+ qS{v x ~ l6D/3)5{v y - WD/3). 



(14) 



Here, the parameter D is chosen to match the initial desired temperature (chosen the same as in the Gaussian initial 
condition). All the quantities are measured every 250 collisions, so that each four consecutive points in figure [1] 
correspond to a time span of one collision per particle. It can be seen that, after some transient, memory of the initial 
condition is forgotten, so that the stationary distribution (/3 = 1) is reached following a universal route. In figure [TJ 
those data points associated to the Gaussian initial distribution approach the scaling curve from above (circles) while 
those for the initial asymmetric case (|14p approach the scaling curve from below (squares). A very similar behavior 
can be seen in Fig [2] for a = 0.8 and four different initial temperatures, again such that To 3> T s , which ensures 
that (3 < 1 (we have also probed the regime /3 > 1 obtained with Tq <C T s , where similar conclusions hold, see e.g. 
Fig. [?] below). We have started either with a Maxwellian distribution, as above, or with a distribution in which all 
the velocities have the same probability density in a square centered on v = (referred to as the "flat" case). Wc 
emphasize that the initial transient is fast: memory of the initial condition is lost after at most 3 or 4 collisions per 
particle, a phenomenon that cannot be appreciated from the figures. 
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Figure 1. Coefficients 0,2 and 03 from Monte Carlo, for a = 0.95. We start with a Maxwellian (circles) or asymmetric (squares) 
distribution, and three initial temperatures: Tq/T s — 22,33,44. Note the vertical scale. 



Borrowing ideas from the extended self-similarity technique [34[, we put to the test the possibility of an enhanced 
universality, by plotting 03 as a function of 02, see Fig. [3] In doing so, it appears that the universal part of the ai 
versus /3 curve seen in Fig. [I] is not enhanced by the reparamctrization 03(02): different initial conditions do not lead 
to a data collapse, beyond the interval —0.22 < a-2 < —0.17 that was already evidenced in Fig. [1] However, a given 
functional form (say Maxwellian) leads to a unique path in the 03-02 plane, which is already a non trivial point, and 
furthermore, a plot like Fig. [3] leads to a significantly reduced scatter of points that Fig. [TJ It is therefore more 
amenable to chart out the universal regime sought for. In these figures, the first measure reported after the dynamics 
has acted on the initial conditions is for a time of 0.25 collisions per particle for the Maxwellian distribution and of 
around 3 collisions per particle for the flat and asymmetric distributions. The present results establish numerically 
the existence of a universal non-trivial scaling regime, for which we now seek analytical characterization. 
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Figure 2. Same as Fig[T] but for a more dissipative system with restitution coefficient a — 0.8. Here, the initial condition is 
either Gaussian (circles, approaching the master curve from above), or flat (squares, see text, approaching the master curve 
from below). In both cases, the initial temperature is To/T a — 7.3, 21.8, 36.4 and 51. 




Figure 3. Same data as in Fig. [T] where the coefficient 03 is shown as a function of 02- In addition to Maxwellian and 
asymmetric, flat initial conditions are also shown. The arrow indicates the steady state values. 

C. Analytical approach 



For the sake of analytical progress, it is convenient to change variables in the Boltzmann equation ([2]), from the set 
{t, v}, to {/3,c}. In these variables, the scaled distribution function fulfills 



[/x(/3) - /z(l)/3 3 ] i — ■ [cix(ci , (3)} + fex(ci , (3) 



dci 



9/3' 



dc 2 f (c 1 ,c 2 )x(ci, ( 8)x(c2, ( 8) + ^(l)/? 3 iLx( 



where 



1 

2d 



dc\' 



dci / dc 2 (^ + ^)To(ci,c 2 )x(ci,/3)x(c2,/3). 



(15) 



(16) 



We note that here and in contrast with Eq. (|T0|). the equation is fully consistent as it appears as a change of variables 
where j3 simply plays the role of time. Nevertheless, proving that for any "reasonable" initial condition, the system 
forgets the initial condition and reaches a universal state is a formidable task. For this reason we limit ourselves to 
the simplified problem of deriving an approximate expression for this distribution function. As in the stationary state, 
the distribution will be worked out in the first Sonine approximation 



X (c^)~XM(c)[l + a 2 (P)S 2 (c 2 )}, 



(17) 
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where the kurtosis, a 2 , has been defined in (|12[) and, by definition, we have 

J dc X (c, 0) = 1, j dcc X (c, 0) = O, J dcc 2 X (c, P) 



(18) 



In expansion (|17l) . we neglect contributions in 03 and higher order. This is justified as long as the inelasticity is not 
too strong, and is backed up here by the fact that |a.3 1 < 02, as can be seen in Figs. [T]and[2] Of course, inclusion of 
higher order terms in (jTTJ) would improve the accuracy of the subsequent calculation. 

Inserting (fT7|) into the Boltzmann equation (|T5l) . taking the fourth velocity moment while neglecting nonlinear 
terms in 02, we obtain the following evolution equation for the cumulant 02 



(1 - B - /3 3 )a 2 (/3) + Ba s 2 



where the parameter B depends on dissipation and space dimension according to 

73 + 8d(7 - 3a) + 15a[2a(l - a) - 7] 



B 



16(1 - a) (3 + 2d + 2a 2 ) + a|[85 + d(30a - 62) + 3a(10a(l - a) - 39)] ' 



(19) 



(20) 



Eq. (TT9")) is an inhomogeneous linear differential equation that can be integrated. It exhibits two singular points at 
/3 = 0, j9 = 1, and we start with the interval [0, 1]. In this case the general solution of the associated homogeneous 
equation reads 



ag(J3)=K 



(1-/3 3 )S J 



(21) 



and a particular solution can be obtained by variations of parameters. The general solution will then be the sum of 
these two contributions. The ensuing a 2 depends on the initial conditions through K, and since our purpose here 
is to extract the universal behaviour of ai as a function of /3, we note that the contribution (|2ip fades rapidly as 
B approaches unity (as (1 — /3 3 ) 45 / 3 where B can be large ; note that it diverges in the elastic limit a — > 1). The 
universal behaviour is consequently encoded is the particular solution, and we finally have 



1 -P 3 
B-l 



1F1 



1 4B - 1 „ 
-- 1- • B 3 

3' ' 3 ' P 



< p < 1. 



where 2-P1 is the hyper-geometric function [32(. This expression is well behaved in all the interval [0, 1] 
analysis can be performed for /3 > 1. Following similar lines, we identify the universal solution to be 



(22) 



An analogous 



a 2 (/3) 



ABa% 



7/3 3 (1-1//3 3 )t 



4B 
1" ; 



10 



1 

w 



(23) 



Clearly, the same technical procedure can be applied to the higher order cumulants. For the sake of simplicity, we 
restrict to the function «2 , that we wish to test against simulation results. 

In order to compare the above theoretical predictions to the simulation data, attention should be payed to the 
fact that analytical computation of velocity moments or cumulants is plagued by non-linear effects that have been 
discussed in the literature pjl, Hlj]. This results in some error in the calculation of the steady state value a 2 , and we 
can also expect B to suffer from a similar inaccuracy, that may be of the order of 10 or 20%. To circumvent this 
(somewhat minor) drawback, we take a| appearing in Eqs. (|22|) and (|23| from the Monte Carlo simulations, and we 
adjust B to match the measured function 02 (/3). This procedure provides us with Fig. |4j where the agreement between 
the functional forms (|2"2"|) and (|2"3"|) with Monte Carlo is excellent. It is of course important to check a posteriori that 
a| and B thereby obtained are close enough to our predictions. The precise values are reported below: 
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(24) 




Figure 4. Kurtosis as a function of /3 (reduced inverse typical velocity), for a system with a = 0.95 (left) and a = 0.80 (right). 
The points are simulation results starting from a Maxwellian distribution and the solid line is the theoretical prediction, given 
by Eq. (|22[) for /3 < 1 and Eq. (|23|l for /3 > 1, see text. The two values of the initial temperature mentioned for each graph 
are used to generate separately the /3 < 1 branch (associated to large To) and the f3 > 1 branch (obtained for small To). The 
steady state values correspond to j3 = 1. 



III. CONCLUSION AND PERSPECTIVES 

To summarize, we have studied the dynamics of a system of inelastic hard spherical grains, heated homogeneously 
by a stochastic thermostat. We have restricted the analysis to low density systems, amenable to a Boltzmann equation 
treatment. We have found that generically, after a kinetic transient, the system evolves into a scaling solution that 
no longer depends on initial conditions, before the steady state is finally reached. The relevant scaling form is not of 
the single parameter family as is the case in the homogeneous cooling state, but involves 2 parameters. The velocity 
distribution function was calculated in the so-called first Sonine approximation, which provides a very good agreement 
with the Monte Carlo simulations. 

At this point several questions arise: What is the counterpart of the universal state brought to the fore at two- 
particle level, or even iV-particlc ? What is the effect of density, and of a change in the driving mechanism ? Do 
the hydrodynamic relations (worked out say at Chapman Enskog level) depend on the structure of this state? The 
complete answer to all these interrogations requires further studies, but we expect that similar scaling forms should 
occur at higher densities, and with different thermostats as long as a steady state can be reached. In this respect, 
the Gaussian thermostat is presumably singular, since it can be mapped onto the free cooling case. The questions 
pertaining to hydrodynamics are more subtle; How the universal /3-scaling behaviour discussed in the present work 
impinges, as a reference state, on transport properties, should be explored. 
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